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Abstract 



Atmospheric aerosols can cause serious damage to human health and life expectancy. 
Using the radiances observed by NASA's Multi-angle Imaging SpectroRadiometer (MISR), 
the current MISR operational algorithm retrieves Aerosol Optical Depth (AOD) at a spatial 
resolution of 17.6x17.6 km 2 . A systematic study of aerosols and their impact on public 
health, especially in highly-populated urban areas, requires a finer-resolution estimate of 
the spatial distribution of AOD values. 
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We embed MISR's operational weighted least squares criterion and its forward simu- 
lations for AOD retrieval in a likelihood framework and further expand it into a Bayesian 
hierarchical model to adapt to a finer spatial scale of 4.4x4.4 km 2 . To take advantage of 
AOD's spatial smoothness, our method borrows strength from data at neighboring pix- 
els by postulating a Gaussian Markov Random Field prior for AOD. Our model considers 
both AOD and aerosol mixing vectors as continuous variables. The inference of AOD 
and mixing vectors is carried out using Metropolis- within-Gibbs sampling methods. Re- 
trieval uncertainties are quantified by posterior variabilities. We also implement a parallel 
MCMC algorithm to reduce computational cost. We assess our retrievals performance 
using ground-based measurements from the AErosol RObotic NETwork (AERONET), a 
hand-held sunphotometer and satellite images from Google Earth. 

Based on case studies in the greater Beijing area, China, we show that a 4.4 km resolu- 
tion can improve the accuracy and coverage of remotely-sensed aerosol retrievals, as well 
as our understanding of the spatial and seasonal behaviors of aerosols. This improvement is 
particularly important during high- AOD events, which often indicate severe air pollution. 

Keywords: Hierarchical Bayesian model, MCMC, spatial dependence, fine retrieval res- 
olution, remote sensing 



1 Motivation 

Atmospheric aerosols, a complex mixture of solid particles and liquid droplets in the air, 
can significantly affect human health and life expectancy. Studies by Monleau et al. and 
Pope et al. show that when inhaled, aerosols can penetrate cell membranes, then mi- 
grate and seriously damage human respiratory and cardiovascular systems [22] and the 
brain [21]. Short-term impacts include irritation to eyes, nose and throat; upper respiratory 



infections including pneumonia and bronchitis; and stroke or death from cardiovascular 
causes. The "Great Smog" in London in 1952, caused by burning coal for heating, left 
approximately 12,000 dead [1]. Continual exposure to hazardous aerosols can aggravate or 
complicate medical conditions in the elderly [4]. Aerosols from silica and diesel can cause 
serious diseases including silicosis and black lung. In addition, aerosols with an aerody- 
namic diameter less than 2.5 /urn, such as black carbon, can severely reduce ground-level 
visibility. Profiling the spatial distribution of aerosols at a fine resolution is critical for air 
quality and public health studies, especially in urban areas with complex anthropogenic 
aerosol sources, such as vehicles, power plants, and factories that burn fossil fuels. 

There are two approaches to measure the spatial distribution of aerosols: through ground- 
based measurements and remote-sensed radiance imageries. Both adopt the same quantifi- 
cation of aerosols, the spectral Aerosol Optical Depth (AOD), the negative logarithm of the 
fraction of radiation not scattered or absorbed on a path 1 . AOD at different spectral bands 
can be viewed as known functions of AOD at the green band using the Angstrom expo- 
nent [20]. For notational simplicity, we refer to AOD throughout this paper as that at the 
green band. With either approach, the spatial and temporal variabilities of aerosols require 
continual observations and computationally efficient analyses. 

The AErosol RObotic NETwork (AERONET) [13] provides a data archive of local 
AOD values using a network of automatic sun photometers (Figure 1, left panel) located at 
more than 400 stations on the Earth's surface. It measures AOD from every half hour to ev- 
ery two hours, with AOD uncertainties < ±0.01 at wavelengths > 440 nm [16]. AERONET 
measurements are widely accepted as the gold standard to validate AOD estimates from 
other data sources. The sparse and heterogeneous locations of AERONET stations, how- 
ever, make it difficult to directly use their measurements to study AOD's spatial behaviors. 

Remote-sensing radiometers offer a better spatial coverage by retrieving AOD from 
radiance imageries over the Earth's entire surface: the Multi-angle Imaging SpectroRa- 



'For example, an AOD value of 2.5 corresponds to 92% of radiation scattered or absorbed. 




Figure 1: AERONET sun photometer at Avignon, France (left) and MISR cameras (right). 



diometer (MISR) aboard the NASA Earth Observing System Terra satellite (Figure 1, right 
panel). MISR views the day lit Earth atmosphere almost simultaneously at nine angles along 
its track. It provides a four- spectral imagery 2 at a l.lxl.l km 2 resolution for all spectral 
bands, but a 275x275 m 2 resolution for the red. 

To quantitatively represent aerosol mixtures, aerosol particles are characterized and cat- 
egorized according to their properties such as radius and single scattering albedo (SSA) 3 . 
A category is referred to as a component aerosol. Then an aerosol mixture is identified by 
its composition: a collection of M component aerosols and their mixing vector relative to 
the M components. Elements of the M-dimensional mixing vector add up to 1, indicating 
mixing percentages of the M components. To simplify remote- sensing retrieval, MISR op- 
erational algorithm considers only 21 component aerosols and 74 pre-fixed compositions 4 . 

Based on these 74 compositions and a discrete grid of AOD values, forward radiative 
transfer (RT) calculations simulate radiances in 36 channels (9 viewing angles x 4 spec- 
tral bands). The results form the Simulated MISR Ancillary Radiative Transfer (SMART) 
Dataset. The MISR operational aerosol retrieval algorithm compares SMART data to the 



2 The four spectral bands are blue, green, red, and near-infrared. More details of MISR observations are 
included in Appendix 7.1. 

3 SSA is defined as the ratio of scattered radiation to total extinct radiation (scattered and absorbed). 

4 Examples of MISR's 21 component aerosols and 74 compositions are displayed in Table 1 and Table 2 
in Appendix 7.2. The non-zero elements of these composition vectors are no more than three. 



MISR observations using a weighted least squares criterion to determine whether they pro- 
vide good fits to the data, i.e., a "successful" retrieval 5 . 

MISR's unique design of multiple viewing angles provides an enhanced sensitivity to 
aerosol scattering and cloud reflective effects [8], giving MISR a significant advantage over 
other remote- sensing instruments. Validated by AERONET measurements, MISR field 
instruments, and airplane campaigns [5], MISR's retrievals have shown to be informative 
in characterizing aerosols. Previous studies include those on wildfire smoke [18], mineral 
dusts [19], and climate changing aerosols [24]. 

MISR's ability to capture aerosol-related information make it well suited to assist stud- 
ies on aerosols' impact on public health. MISR operational retrievals at the 17.6 km reso- 
lution, however, are too coarse to determine the relationship between aerosols and human 
health. The heterogeneity of urban aerosols within an area of 17.6x17.6 km 2 makes a finer 
resolution desirable. For example, San Francisco is represented by less than half of a pixel 
on MISR output imagery. Yet the residents of San Francisco are exposed to varying level 
of air pollution. One also needs to add more varieties beyond the 74 pre-fixed compositions 
to capture the growing heterogeneity of urban aerosols. 

Finer-resolution retrievals with greater varieties of aerosol compositions lead to a larger 
number of parameters for us to estimate. This retrieval is possible if we take advantage of 
AOD's spatial smoothness and reduce the collection of 21 component aerosols to a smaller 
subset, say four. We chose to take a MISR Block 6 as a data unit to balance the coverage of 
a greater metropolitan area and computational cost. 

Our study proposes a Bayesian hierarchical model to retrieve AOD values and mix- 
ing vectors given a collection of four component aerosols at a resolution of 4.4x4.4 km 2 . 
The retrievals are based on MISR observations. The four component aerosols are chosen 
according to current knowledge of aerosol conditions for the study region. We adopt a 



5 More details on MISR operational aerosol retrieval algorithm are included in Appendix 7.2. 

6 MISR observes the Earth's surface in 233 swaths; each swath contains 180 560x140 km 2 MISR Blocks. 



likelihood framework based on MISR's weighted least squares and construct a Bayesian 
hierarchy to incorporate AOD's spatial smoothness into the model. In particular, we intro- 
duce a Gaussian Markov Random Field (GMRF) prior to model AOD's spatial dependence 
structure. The movement and dispersion of air particles in the atmosphere justify this spatial 
smoothness of AOD from a physical viewpoint. To flexibly describe various aerosol con- 
ditions, our model regards AOD values and mixing vectors as continuous parameters. This 
expands the variety of possible compositions beyond the 74 MISR-designated choices. We 
show how this enriched variety is necessary for retrieving heterogeneous urban aerosols. 
We have not introduced an explicit spatial dependence structure for the mixing vectors of 
pixels in one MISR Block. Instead, the spatial smoothness of mixing vectors is implicitly 
enforced by the smoothness of AOD values through the hierarchy. We introduce conjugate 
hyperpriors to determine the parameters in the priors. 

The posterior inference of AOD and mixing vectors is carried out using Markov Chain 
Monte Carlo (MCMC) sampling methods, particularly Metropolis- within-Gibbs. The sam- 
pling methods also allow us to quantify the retrieval uncertainties by posterior variability. 
The algorithm, however, is computationally intense. We develop a parallel MCMC algo- 
rithm by breaking the MISR block into smaller patches to enable parallel samplings while 
maintaining the overall smoothness using summary statistics. We show that retrievals from 
the two algorithms are consistent, with an increase in computational speed for the latter. 

To assess the performance of our methods, we apply them to retrieve AOD and mixing 
vectors for the greater Beijing area in China. Our retrievals are tested against ground-based 
measurements from two AERONET stations in the area. They show improvement on ac- 
curacy and retrieval coverage, especially in high- AOD events. We also include geographic 
information from Google Earth to qualitatively validate our retrievals. 

The rest of the paper is organized as follows: Section 2 provides the rationale and de- 
tails of our Bayesian model for retrieving AOD values and mixing vectors, while Section 3 



details our MCMC algorithms. Section 4 contains case studies for model validation and in- 
terpretation, comparing our results with MISR's retrievals and AERONET measurements. 
In particular, Section 4.4 uses Google Earth to investigate locations with extreme AOD 
values, relating the retrievals to local geophysical conditions and levels of anthropogenic 
activities. Section 4.5 illustrates the necessity to include a richer variety of compositions. 
Section 5 summarizes the results and suggests directions for future research. 

2 Hierarchical Bayesian Model 

Our study objective is to establish a more detailed data-driven description of the relation- 
ships among radiances, AOD, and aerosol compositions to assist aerosol-related health 
studies. The MISR operational retrieval algorithm provides this information by comparing 
the observed and RT calculated radiances, but it is limited within the 74 pre-fixed compo- 
sitions and a discrete grid of AOD values. We propose to understand a greater variety of 
optical behaviors of aerosols by considering AOD values and mixing vectors as continu- 
ous variables, given a fixed collection of four component aerosols. For the greater Beijing 
area, the collection includes spherical non-absorbing aerosols without sulfate, spherical 
non-absorbing aerosols with sulfate, spherical absorbing aerosols, and grains (dust). 

Each MISR Block contains 256 pixels (8 rows x 32 columns) at the 17.6 km resolu- 
tion in MISR retrievals. The number of pixels rises to 4,096 (32 rows x 128 columns) as 
we retrieve at 4.4 km resolution, presenting a more complex problem with approximately 
16,384 parameters to estimate for each of 4 1,940 MISR Blocks. On the other hand, air par- 
ticles interact in the atmosphere within a certain range; consequently, they affect aerosol 
conditions in near neighborhoods [2] [26]. This suggests a stronger spatial dependence 
among adjacent pixels at a finer scale. When modeling at a fine resolution, therefore, it 
is beneficial and necessary to borrow strength from AOD's spatial smoothness to reduce 
model complexity. In particular, we construct a hierarchical Bayesian model with a built-in 



spatial dependence structure using a GMRF prior for AOD values. 

In Section 2.1, we introduce the notations and define the likelihood function, based on 
MISR's weighted least squares criterion. Section 2.2 introduces our Bayesian hierarchy. 
Unless specified otherwise, our Bayesian retrievals are at a 4.4 km resolution. 

2.1 Defining the Likelihood Function 

Let p = 1 , . . . , P denote the P = 4, 096 pixels on a two-dimensional lattice in a MISR Block 
at 4.4 km. For pixel p, MISR observes the top-of-atmosphere radiances L p = (L\ p , . . . , L Cp ) 
at MISR's C = 36 channels. Every channel c is weighted by -K, c = 1, . . . ,36, to ac- 
count for different measurement errors. For pixel p, let t p e R denote its AOD value, 
and 6 P = (6 p i, . . . , 6 pM ) e R M the mixing vector of the M component aerosols involved. 
Zm=i°pm = 1. Let L = {L u ...,Lp}, t = {ti,...,t p }, and 6 = {0i,...,0 P }. For each 
pixel p independently, the RT calculated radiances L RT = (L RT , . . . , L RT ) can be viewed as 
functions of (r p , 6 P ), relative to the collection of component aerosols involved. The MISR 
operational retrieval algorithm exhaustively searches over combinations of pre-fixed AOD 
values and 74 compositions to match L RT to the observed, based on the following weighted 
least squares [6]: 

2 _ y ( L cp ~ L c ( T /» A/>)) m 

The combinations of AOD and compositions whose xp satisfy a pre-established threshold 
are considered successful and the averages of all successful retrievals are the outputs. 

In this paper, we turn the weighted least squares in (1) into a likelihood function to 
assess how probable are the observations L given r and 6. We quantify this probability 
through a Gaussian density based on the Euclidean distance between L and L RT . Aerosol 
retrievals constrained by (1) can be formalized as inference of (r, 6) based on the following 



likelihood function, 

rtHr.W^pj-EE ^'-y'" 2 }. O 

In the MISR operational retrieval, {cr^}^ =1 are explicitly assigned as the minimum of 0.04 
and L c = (£p=i L cp )/P (Appendix 7.2). In our model, {o"^}^ =1 are considered unknown 
parameters and are to be estimated together with r and 6. We set the number of component 
aerosols M = A and case studies show the sufficiency of this choice in our cases. 

We adopt a Bayesian approach since it easily incorporates smoothness in r through their 
prior distributions. It also provides uncertainty measures through posterior variances. In the 
following section, we describe our hierarchical model by assigning reasonable priors to the 
unobserved variables and building conditional relationships within the Bayesian hierarchy. 

2.2 Construction of the Priors and Conditional Probabilities 

For fixed atmospheric pressures, humidity and wind levels, the top-of-atmosphere radiances 
L are mainly determined by r and 6. We use the hierarchy's first level to indicate the 
dependence of L on r and 6. Prior distributions for r and 6 are further assumed to be 
independent, i.e., p(r, 6) = p(T)p(6) to simplify computation. The posterior distribution of 
r will influence that of 6 by interacting within the hierarchy, and vice versa. 

Priors for r and 6 are postulated to capture their smoothness calibrated by hyperparameters 
k and a respectively. A GMRF prior is adopted to describe r's spatial dependence. Though 
there is no explicit spatial structure for 6, the interactions between r and 6 through the 
hierarchy implicitly enforce a spatial dependence for 6. Figure 2 provides a graphical rep- 
resentation of the hierarchy. It characterizes the dependence of L on parameters r and 6 
and the independence among the hyperparameters. The inference of parameters and hyper- 
parameters, using MCMC sampling methods, is discussed in Section 3. 




Figure 2: Graphical representation of the hierarchical Bayesian model. Green square: ob- 
servations; red circle: model parameters; blue circle: model hyperparameters. 

2.2.1 Prior Beliefs about AOD's Spatial Dependence 

We model the spatial dependence structure of AOD values using an intrinsic GMRF prior 
of first order [23]. Define k as the homogenous scaler precision and use ~ to indicate spatial 
adjacency. The following prior is invariant to perturbation by the same constant to r of all 
pixels [3], 

P (t\k) oc iC^ exp | - ^ Yj ( V - t p) 2 [ ■ ( 3 ) 

This allows us to model the smoothness of r, regardless of its unknown overall level, by 
penalizing sharp changes in r among adjacent pixels. 

The above GMRF prior is governed by one precision parameter k. The larger k is, the 
smoother the region's AOD values are, with 4= as AOD's precision. For some regions, 
however, a more complicated GMRF prior might be necessary. For example, a continu- 
ing pattern of wind directions might make it useful to distinguish a upwind pixel from a 
downwind pixel. This paper works with a homogenous precision k for r's GMRF prior. 

To estimate k, we assign k a hyperprior. Since p(t\k) belongs to an exponential fam- 
ily, we choose its natural conjugate as p(k), a noninformative Gamma distribution. The 
posterior is also a Gamma distribution, 

P (k\t) <x k (P ' 1)I2 ' 1 exp{-^ Yj (V - t p^- W 

p':p'~p 
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2.2.2 Prior Specification for Aerosol Compositions 

Prior information on aerosol compositions is incorporated in the model through the choice 
of the M = 4 component aerosols involved, based on geophysical knowledge of the study 
region. To model the mixing vectors 9 for these M component aerosols, we adopt an M- 
dimensional Dirichlet prior with parameter a = (a\,..., ckm)- By adjusting the magnitude 
of a, we can control the sparsity of the mixings of component aerosols. Conditioning on 
a, the mixing vectors 9 P for each pixel p are considered to be independent, 

P P rvY M \ 

p(9\a) = f] p(9 p \a) = f] \1^ X ^ ■ ■ ■ ^ ■ (5) 

p=\ p =l llm=l l \ a m) 

In general cases with no information on the mixing's sparsity, we assign a a hyperprior 
to estimate it. Since (5) belongs to an exponential family, we take its conjugate as p(a), 

/ M m m \ 

p(a) oc exp ^(a m - l)y ilm - rj(^ log T(a m ) - log T(^ a m )) . 

\m=\ m=\ m=\ ) 

7] is chosen as and {yo, m }ff =1 as 1 to obtain an exponential hyperprior. The posterior is, 

MP M M 

p(a\0) oc exp{^(a m - 1)(^ log 6 pm + l)-Px (^ log T(a m ) - log T{^ a m ))}. 

m=\ p=\ m=l m=l 

2.2.3 Hyperpriors for cr 

The likelihood function for cr = (en,. . . ,cr c ), p(L\cr,T,9), follows a normal distribu- 
tion with known mean and unknown variance. We adopt a conjugate noninformative 



scaled inverse-^ prior for cr to model the channel weights {tM. Thus, p(cr 2 .) oc cr c 



2 



(c = 1, . . . , C), and the conditional posterior also follows the scaled inverse-^ 2 distribution, 



P(o-1\t,0,L)cc(o- 2 c ) ( 2 +1) exp<^ '- 



2al 
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3 MCMC Retrieval Algorithms 

Based on the hierarchical Bayesian model previously developed, this section derives marginal 
posterior distributions of AOD values r and mixing vectors 6 and devises two MCMC algo- 
rithms to sample from the posteriors. Using MISR observed radiances as input, the sampled 
posterior means are used as our retrieval outputs. 

3.1 Posterior Distributions of AOD and compositions 

The full Bayesian model discussed above can be summarized as follows: 



l p \t p , e p ~ n(l rt (t p , e p ), (t\ p = i,...,p, 



t\k ~ GMRF{k), 
9\a ~ Dirichlet{a), 

(j ~ scaled inverse ~x 2 ( v o)> 

k ~ Gamma(aQ,{3 ), 

M 

p(a) ~ Expi^il - a m )). 



m=\ 



With no additional information on the hyperparameters, v , a , and /3 are chosen to be 
for convenience and later shown to be robust. The marginal posterior of r is, 



p(T\9, k, a, L) oc exp Uk £ <J, - r p f - £ £ ^ cp - L?(t p ,9 p )? ^ (g) 

I p'-p'~p 



2(r ' 



The marginal posterior distribution of 6 can be expressed as, 

p(0\r, a, cr, L) ^ exp J ^ £(<*„ - 1) log 9 pm - J] J] ^EEZR^ElM I . (7 ) 

\p=l m=\ c=l p=\ c ) 
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Both posteriors contain an RT-based L RT , can be obtained at necessary values through inter- 
polations from MISR SMART Dataset [7], using t p and 6 P as input (Appendix 7.4). These 
non-closed-form posteriors, however, are difficult to directly sample from. A Metropolis- 
within-Gibbs sampler is used when RT-based L RT are involved. 

3.2 Metropolis-within-Gibbs Sampling from the Posterior Distributions 

The Gibbs sampler [11] is a numerical technique to sample from a joint posterior distribu- 
tion, p(T,9,cr,K,a\L) in our case. We sample for t p and 6 P using a Metropolis-Hastings 
(M-H) sampler, for each pixel p on the MISR Block column-by-column and pixel-by-pixel. 
The following proposal distribution is used in M-H for t p , 



n p K 1 -r-, 

P(t p \t. p ) oc exp — — (t p 2j V- 

V p p':p'~p 



where n p is the number of adjacent pixels to pixel p. A Dirichlet proposal distribution with 
parameter a is used for 6 p . Denote vector (t p , . . . , t p >) by r p:p > and similarly for 6 and a. 
Given initializations (r (0) , 6 <0> , cr (0) , k (0) , or (0) ), the Gibbs sampler proceeds as follows: 
Metropolis-within-Gibbs Algorithm (M-w-G) 



At step t, iterate the following process: 
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The above cycle generates a realization of a Markov chain, which gives approximate 
samples from the marginal posteriors after a successful burn-in process [11]. Appendix 7.5 
contains details on the stability of sampled posterior mean over different initializations. 

Whenever we used the Metropolis-Hastings algorithm, we checked that the acceptance 
rate was roughly between 25% and 50% to secure adequate mixing of posterior samples 
[9]. We also used a variety of diagnostic tools to ensure the convergence of the Markov 
chains, including the potential scale reduction R [10] on several parallel Markov chains 
with different initializations 7 . We ran the chains until R is less than 1.1 or 1.2, using R of 
the logarithm of the posterior distribution as a benchmark. In addition, a geometric decay 
of the autocorrelation as a function of the lag also validated the mixing of our chains. 

In the end, we estimate the posterior mean based on the samples and consider it the 
point estimates for r and 6. While the MCMC algorithm enables us to handle a hierarchy 
whose complexity precludes fitting by analytical methods, its computation intensity limits 
its online use. Next, we propose a parallel MCMC algorithm to reduce computational cost. 

3.3 A Parallel MCMC Algorithm 

Many MCMC sampling algorithms for spatial data suffer from high computational cost 
caused by the large dimensionality of spatial data. On a 4.4 km scale, our MCMC algo- 
rithm simulates samples for more than 16,000 variables 8 for one MISR Block. The large 
computation cost is exacerbated by the non-closed form of the posterior distributions. It is 
crucial, therefore, to develop an efficient and reasonably accurate algorithm. 

We devise a parallel MCMC algorithm to improve the computational efficiency. In 
this algorithm, each MISR block is divided into 2x8 patches of equal size with at least 
four overlapping columns and rows between adjacent patches. The M-w-G sampler is 
applied to each patch independently to generate samples of r and 6 for these patches. This 



7 For example, cases are inspected with initial values of k as 10, 100, 200, 800. 

8 Excluding cloudy pixels can sometimes reduce the total dimensions to around 5,000 for one MISR Block. 
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independent sampling on different patches can therefore benefit from parallel computing. 

Information on AOD's spatial dependence structure is to be maintained across the en- 
tire MISR Block to achieve an accurate estimate of its spatial smoothness level. On that 
account, we let the parallel MCMC algorithm periodically exchange spatial smoothness 
information among all patches on the entire MISR Block. Given /c's conditional posterior, 



P (k\t) <x k {P - 1)I2 - 1 exp{-^ J](T P - v) 2 }, 



2 
p~p 

and summary statistic, T K = Yj P ~ P '( t p ~ t p >) 2 , it follows that p(k\t) = p(k\T k ). Hence T K is a 
sufficient statistic that holds all information on calibration k for AOD's spatial smoothness 
across the entire MISR block. We now describe the parallel MCMC algorithm in detail: 
Parallel MCMC Algorithm 

1. Obtain a standard MISR block of 32x128 pixels. Divide the block into 2x8 patches, each 
of 20x20 pixels, with at least 4 overlapping columns/rows between adjacent patches. 

2. At step t, iterate the following process: 

- Use M-w-G algorithm to sample r ~ p(r\0, k, <x, L), ~ p(0\T, a, <x, L), cr ~ 
/?(cr|r®), k ~ p{k\T^), a ~ p(a\T^) within each patch in parallel for 50 iterations. 

- Average the samples of the overlapping pixels between any two adjacent patches. 

- Calculate summary statistics using current samples, 

t('+1) _ \-P (j jRTr 



* = 2p=i(^/> - L c 0> W. c = i, . . . , c, 

jp~p'\tp ~ T p') 



T i,+l) = y ,(t -t ,) 2 



7t +1) = Zp=ilog0 pm ,m=l,...,M. 



Given that the hyperparameters control the level of spatial smoothness of model param- 
eters in all patches, the spatially dependence structure of this entire imagery is preserved, 
while the parallel-patch samplings improve the algorithm's computation efficiency. This 
parallel MCMC sampling scheme can be generalized for other computations on spatial 
data using sufficient statistics to maintain a global spatial structure. The above process can 
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be automated using the Perl programming language. For one MISR block at a resolution 
of 4.4 km, the computational time for the parallel MCMC algorithm is less than one-fifth 
of that for the global MCMC sampling algorithm (Appendix 7.5). 

The global and the parallel MCMC algorithms produce reasonably consistent results. 
Figure 3 shows two scatterplots of AOD retrievals from global and parallel MCMC al- 
gorithms after the same number of iterations. The outputs generally agree, except for a 
small group of pixels that mostly lie on the edges of patches. The spatial smoothness is in- 
terrupted between patches; therefore, the benefits of a stabilizing factor from neighboring 
pixels are lost. This confirms the importance of maintaining an appropriate spatial de- 
pendence structure. Increasing the number of iterations and communications of summary 
statistics, as well as moving the edges of patches, reduces the disagreement. 

The next section evaluates the performance of our Bayesian retrievals' using case stud- 
ies. For non-operational model validations, we apply the global MCMC algorithm to avoid 
inconsistency in number of iterations for different retrievals. 




0.25 0.3 0.35 




Figure 3: Two examples of the two MCMC algorithms, x: AOD retrievals from the global 
MCMC sampling algorithm; y: AOD retrievals from the parallel MCMC algorithm. 
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4 Assessing MCMC Algorithm Performance: Case Studies on Aerosol 
Retrievals for the Greater Beijing Area, China 

In this section, comparison studies to MISR outputs are presented for the greater Beijing 
area and their differences discussed. We validate our AOD retrievals using AERONET 
measurements and satellite images from Google Earth. Through case studies, we demon- 
strate the importance of high-resolution retrieval to gain more information on aerosols, 
and a greater variety of compositions to improve retrieval accuracy and coverage. All 
retrievals in this section are based on MISR observations over the greater Beijing area (lat- 
itude 38.95N~40.15N; longitude: 115.57E~119.50E). 

4.1 Historical and Present Background of Air Pollution in Beijing 

There are two main reasons for high AOD values in Beijing. First, anthropogenic activi- 
ties produce various air pollutants such as waste gas from over two million automobiles, 
industrial emissions, and dust from over 5,000 construction sites. Coal-dominated energy 
supplies burn off over 28 million tons of coal per year. Due to the presence of multiple air 
pollutants, Beijing bears high AOD and heterogenous aerosol compositions. Second, geo- 
graphical characteristics of the region also contribute to the buildup of aerosols. Mountains 
and hills surround the greater Beijing area on the east, west, and north. Aerosols accumu- 
late until removed by rain or dispersed by wind into the plains on the south. Wind in any 
other direction is blocked and venting of aerosols highly restricted. 

The anthropogenic behaviors and seasonal weather patterns lead to a yearly trend of 
air quality in Beijing. In winter, starting from mid-November, the use of central heating 
systems bring large amounts of sulfur dioxide, nitrogen dioxide, and other types of incom- 
pletely burned waste into the air . Enormous boilers of the heating systems generate masses 
of water vapor that help to increase the absorption of particles in the air. In December, the 
cold airflow from Siberia generates northwestern winds in Beijing that clear the accumu- 
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lated pollutants, leading to the lowest aerosol loadings in January. The loadings vary due 
to the discontinuous cold airflow. In spring, wind from Siberia and Loess Plateau of North- 
western China brings sandstorms. Traffic and industrial activities revive as temperatures 
warm up. From June to August, aerosol loadings reach the highest of the year [25]. From 
September to December, precipitation disperses and dissolves air pollutants [14]. 

Significant measures have been introduced to control air pollution in Beijing [15], but 
rapid economic and population growth, as well as the continuous expansion of the trans- 
portation system, have hindered their success. On March 20, 2010, Beijing experienced 
its biggest sandstorm since 2006. The storm originated from Inner Mongolia and haunted 
Beijing for over 5 hours. Characterizing AOD's spatial distribution under such extreme air 
conditions is especially important to understand its influence on public health. 

4.2 Comparison with MISR Retrievals 

To retrieve AOD, we average all possible AOD values weighted by the posterior distribu- 
tion, while MISR averages successful retrievals with equal weights. This suggests consis- 
tency between our retrievals and those of MISR, and possible differences. 

For a first impression of the aerosol conditions in this region on Jan 10, 2009, Figure 4 
panel (a) displays the observed MISR imagery by one of the forward cameras (CF). MISR 
AOD retrievals are plotted at their 17.6 km resolution in panel (b). Our Bayesian AOD 
retrievals plotted in panel (c) at 4.4 km. Figure 4 suggests shared information in MISR 
and our MCMC retrievals, such as the coastline on the right, the overall AOD level, and 
its spatial patterns. This consistency is confirmed by the scatterplots of MISR outputs and 
Bayesian AOD retrievals aggregated to a 17.6 km resolution (Figure 5, left panel). 

The increased diversity in our Bayesian AOD retrievals across the MISR Block (Figure 
5, right panel) indicates that more details are detected as the spatial resolution improves. 
This is expected, since a finer resolution leads to more information observed and piped 
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(a) MISR radiance imagery by the nadir camera on Jan 10, 2009. 




(b) MISR AOD retrievals at 17.6 km resolution. 




(c) Bayesian AOD retrievals using MCMC at 4.4 km resolution. 
Figure 4: AOD estimates from MISR and our Bayesian retrievals. 

into the model. But the reliability of such diversity needs to be further validated by other 
independent sources such as ground-based measurements (Section 3.4). 

The black pixels represent missing retrievals due to MISR instrumental design or mask- 
ing criteria such as cloud screening. Aerosol retrievals are not attempted when clouds are 
detected. MISR averages the l.lxl.l km 2 observations into a pixel of 17.6x17.6 km 2 and 
ignores small clouds, given the cloudless areas are more than -^ of the pixel. Small clouds 
negligible on the 17.6 km scale might be significant on a 4.4 km scale, so we tend to have 
more missing retrievals in some areas intrinsically determined by the observations. How- 
ever, when none of the 74 MISR-designated compositions produce successful retrievals, 
MISR operational algorithm marks it as missing. Our Bayesian retrievals, with a richer 
variety of compositions, eliminate these unnecessary missing retrievals (Section 4.5). 

To objectively evaluate the performance of our MCMC retrieval algorithms, the next 
section conducts a validation study using ground-based measurements. 
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Figure 5: Scatterplots of MISR against MCMC retrievals at an aggregated 17.6 km resolu- 
tion (left, r.m.s. = 0.0295) and a 4.4 km resolution (right, r.m.s. = 0.0309). 

4.3 Model Validation for Bayesian Retrievals by Ground-based Measurements 

Ground-based measurements are collected at AERONET Beijing and AERONET Xianghe 
stations, as well as via a hand-held MICROTOPS II Sunphotometer at several locations in 
urban Beijing. The fixed locations of the AERONET stations and the limited travel range of 
the Sunphotometer' s human operator make it impossible to validate retrievals of all pixels 
on the MISR Block. Instead, we focus on the pixels that contain the AERONET stations 
and our Sunphotometer- visited locations. To match the Bayesian retrievals at 550 nm, 
we first convert AERONET data to a wavelength of 550 nm using AERONET estimates 
of Angstrom exponent and a linear transformation. We then average the measurements 
within a one-hour window when Terra carrying MISR passes over the AERONET stations. 
Jiang, et al. [17] show that a narrower time window better captures the correlations between 
AERONET measurements and MISR retrievals. The frequent cloudy weather in the greater 
Beijing area and its latitude 9 contribute to the difficulty of matching remote-sensed versus 
ground-based data pairs for model validation. 



9 The Beijing city is visited by the Terra satellite every five to nine days. 
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4.3.1 Retrieval Validation at AERONET Beijing Station 

Figure 6 shows a boxplot of our Bayesian AOD retrievals for the pixel which contains the 
AERONET Beijing Station 10 , with estimated uncertainties indicated by the box edges for 
inter-quartile ranges of posteriors and the whiskers drawn to the 5th and 95th percentiles. 
The three retrievals on March 15, April 30, and May 16, 2009 are plotted separately in the 
right panel to keep an appropriate scale for the left panel. 

As long as a pixel is cloudless, our MCMC algorithms provide an AOD retrieval. How- 
ever, the MISR operational retrieval algorithm produces missing values for 24% of the 21 
cases in Figure 6. This results from the limited choices of aerosol compositions in MISR 
retrievals and the increasingly heterogeneous aerosol conditions in Beijing. 

In general, our Bayesian retrievals show improvement in accuracy. Detailed informa- 
tion on aerosols are revealed by the fine-resolution retrievals. In the coarse-resolution re- 
trievals, however, high AOD values are averaged down by its neighbors and low AOD 
values averaged up, resulting in the loss of useful information. The three high AOD values 
in the right panel of Figure 6 indicate Beijing's extreme air conditions, corresponding to 
86%, 71%, and 81% reduced radiation by aerosols. For example, records of news on the 
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Figure 6: Validation of our AOD retrievals by measurements at AERONET Beijing Station. 
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Latitude: 39.97689° North; longitude: 116.38137° East. 
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Figure 7: Left: air conditions in Beijing, China, on March 15, 2009; right: locations of 
Hebei Province and Shenyang. 

Internet indicate that on March 15, 2009, the city was trapped in a sandstorm originated in 
Inner Mongolia (Figure 7). 

4.3.2 Retrieval Validation at AERONET Xianghe Station 

Figure 8 compares the remote- sensed retrievals to AERONET measurements at the pixel 
that contains the AERONET Xianghe station 11 . From December to February, AERONET 
measurements are mostly higher than remote-sensed retrievals, but no pattern afterwards. 

AERONET Xianghe station has the Jingshen Expressway to its north, which is a major 
path connecting two hub cities: Beijing and Shenyang 12 . The northwest wind in winter 
carries car exhaust to the AERONET Xianghe station, possibly leading to high AOD mea- 
surements. Yet for remote- sensing retrievals, the green fields in a larger neighborhood bal- 
ance this factor, which possibly results in a washed-out signal. However, the fine-resolution 
retrievals seem to suffer less from the balancing factors and display a better accuracy. 

We would like to discuss one of the cases where our AOD retrieval is much higher than 
AERONET measurement: the first data point in Figure 8 (December 25, 2008). MISR 
operational retrieval algorithm produced no output for this day. To the east of Xianghe 



"Latitude: 39.75360° North; longitude: 116.96150° East. 

12 The capital and largest city of Liaoning Province in Northeast China. 
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Figure 8: Validation of our AOD retrieval results by AERONET measurements in Xianghe. 

station in Hebei Province lie several major malls for furniture exhibition and manufacture. 
On December 25, 2008, the furniture companies started renovating their exhibition halls. 
The construction could have caused localized aerosol loadings not observed by the Xianghe 
AERONET site 2 km away upwind within one day, but detected by the MISR instrument 
and captured by our retrievals. 

4.4 Validation Using Google Earth 

Model validation by ground measurements is reassuring, except for the limited number of 
matching data pairs. For areas where ground measurements are not available, we investi- 
gate using satellite images from Google Earth. 

Figure 9 shows a projection of our retrievals onto Google Earth. We take advantage 
of the interactive feature of Google Earth to zoom in on the map to investigate the geo- 
graphical and anthropogenic conditions of each region, for example, the existence of heavy 
industries and transportation patterns. In particular, we examine pixels where the regional 
AOD values are not very smooth. Information carried by these pixels are likely to be re- 
vealed only at a fine resolution. We are able to identify possible emission sources around 
several locations with high AOD values as marked in Figure 9: 
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Figure 9: Retrieval results projected on Google Earth. 

A) A hub of several major highways, including Jingshen and Jingtang Highways; 

B) Industrialized areas, with shadows of high chimneys; 

C) Construction site pollution: building a new holiday resort. 

Similarly, our investigation found locations with a low level of AOD values reveals that 
they have more natural features and less traffic and industries, for example: 

D) The Olympic Park; 

E) Mountain Xi, a famous hiking destination; and 

F) Beidaihe, a famous beach resort. 

Google Earth, with its interactive character, can be employed to understand the behav- 
iors of aerosols and identify locations with high or low AOD values. 



4.5 Case Study for Including a Richer Variety of Aerosol Compositions 

This section emphasizes the necessity of expanding MISR-designated 74 aerosol compo- 
sitions. This expansion improves coverage of the MISR operational retrievals and detects 
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more features of aerosol behaviors (e.g., seasonal trend for certain component aerosols). 

MISR AOD (558 nm) 




MCMC AOD (558nm) 




Figure 10: Case study of AOD retrievals on March 15, 2009. 

For an example of the improvement on retrieval coverage, we examine March 15, 2009. 
MISR failed to retrieve AOD for the majority of the Block (Figure 10). Our Bayesian 
retrievals provide better coverage and the retrievals give distinctly high AOD values with 
a clear path of aerosols migrating from west to east and into the ocean. This unusual 
discrepancy leads us to run through the weather records: on that day, the area suffered from 
a sandstorm originated in Inner Mongolia, which later passed into eastern China. For areas 
like Beijing, which experience occasional sandstorms, the limited compositions containing 
grains(dust) among MISR's 74 choices could easily result in a low coverage of MISR 
retrievals. Similar situations might exist for other locations with usual aerosol conditions. 

The retrieved mixing vectors also contains information on the regional aerosol compo- 
sition and can be used to identify pollution type and source. Figure 1 1 shows the retrieved 
mixing vectors on January 10, 2009, on the left, and those on April 14, 2009, on the right. 
Component No. 6 with sulfate tends to dominate the composition in winter due to coal 
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Figure 11: Retrievals of aerosol composition with = (.61,62,63,04) plotted as lengths of 
segments in a column. Left: January 10, 2009; Right: April 14, 2009. 

burning for heating, with No. 19, grains (dust) dominating in spring due to sandstorms. 

Figure 12 shows the mixing percentages of component No. 19 over December 2008 to 
June 2009, at four different locations: the AERONET Beijing station, the AERONET Xi- 
anghe station, location (A) and (F) marked in Figure 9. For AERONET Beijing station, the 
percentage of grains(dust) only rose in the spring, due to the sandstorms, while the con- 
structions around AERONET Xianghe might have raised the percentage earlier in the year. 
Location (A), where major highways intersect, showed a high amount of dust in its aerosol 
compositions through the warm seasons when traffic typically increases. The mixing per- 
centage of No. 19 at Beidaihe Resort moved relatively in consistence with AERONET Bei- 
jing station. We hope to explore this trend in future research. 

In general, by correctly identifying the major pollutants for each season, we can better 
understand the transitions of aerosols and, therefore, take efforts to improve air quality in a 
more specific and to the point manner. For accuracy and coverage, it is necessary to expand 
the MISR-designated 74 compositions to a richer variety. 
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Figure 12: Mixing percentages of component No. 19 from winter to spring. 

5 Discussion 

Aerosols serve as an important factor in climate systems and air quality. Global climate and 
regional wind patterns, as well as anthropogenic activities, shape aerosols into an increas- 
ingly heterogeneous system. Such complexity calls for flexible models with a compatibly 
fine resolution. In this paper, we have presented a Bayesian hierarchical model to retrieve 
AOD and mixing vectors relative to a collection of four component aerosols at a spatial res- 
olution of 4.4 km using MISR observations. The model incorporates a spatial dependence 
structure to gain strength from AOD's smoothness, at the same time allows for a richer 
variety of aerosol mixing vectors. A more detailed AOD spatial profile is provided then 
validated by AERONET and Google Earth; a better retrieval coverage is obtained due to 
the flexible choices of aerosol compositions and improved accuracy. This improvement is 
particularly important during high- AOD events, which often indicate severe air pollution. 

We hope that a more detailed aerosol profile will eventually expand the potential of 
MISR observations to assist urban air quality monitoring and public health studies. From 
the cases, we become more aware of the complexity in modeling aerosol conditions. Ex- 
ternal data such as Google Earth can introduce new access to model validation. 

For future research, we plan to incorporate more prior information into the model, such 
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as wind measurements. If this information is available on an appropriate scale, it may 
improve retrieval performance. Our ultimate goal is to explore the possibility of using our 
AOD retrievals at a fine resolution to help understand the impact of AOD on public health. 
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7 Appendix 

7.1 The MISR Instrument and Observations 

MISR views the daylit Earth atmosphere almost simultaneously at nine angles along its 
track. It provides a four-spectral imagery at a 275x275 m 2 spatial resolution. The radiances 
in all spectral bands except for the red are aggregated onboard to a l.lxl.l km 2 resolution 
to reduce data transmission from Terra. The MISR cameras cover the Earth's surface in 233 
geographically distinct but overlapping paths of approximately 560 km wide, extending 
from the Arctic to Antarctica. Every swath is divided into 180 140-km-long blocks. 
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7.2 The MISR AOD Retrieval Algorithm 

In this section, we describe the MISR Level 2 aerosol retrieval algorithm [6]. Let p = 

Table 1: Examples of single component aerosols, 21 in total. 
No. Name Radius (fxm) SSA (558nm, Green band) 

I spherical_nonabsorbing (sulfate, sea salt, organic) 0.0562 1 

II sphericaLabsorbing (sea salt, organic) 1.2835 0.899 
19 grains_model (dust) 0.7542 0.977 

1, . . . , P denote the P pixels on a two-dimensional spatial lattice, i.e. a MISR block. For 
every pixel p, MISR observes the top-of-atmosphere sun reflectance L p . Each L p is a 
C = 36 dimensional vector, corresponding to MISR observations for pixel p at the 36 
parallel instrument channels: the four spectral bands in each of the nine viewing angles. 
Every channel c is weighted by cr c , c = 1, . . . , 36, to account for different measurement or 
modeling errors. In MISR retrieval algorithm, cr c is assigned as, cr c = 0.05 xmax(L c , 0.04). 
where L c is the average of observed sun reflectances at channel c over all pixels. 

Based on a balance between computation expense and current knowledge on aerosols 
mixings in various geographical locations, MISR retrieval algorithm only considers aerosol 
compositions containing M = 3 or fewer component aerosols, examples shown in Ta- 
ble 2. Our model will extend to allow aerosol compositions containing M = A com- 
ponent aerosols. For example, the four components involved in aerosol retrieval for the 
greater Beijing area include spherical non-absorbing aerosols (without sulfate), spherical 
non-absorbing aerosols (with sulfate), spherical absorbing aerosols, and grains(dust). 

For the pth pixel on the lattice, use t p to denote its AOD value, and 6 P = {6 P \, . . . , 8 p m) 
for the aerosol composition, i.e. values of the associated M-dimensional vector indicating 

Table 2: Examples of aerosol compositions, 74 in total 

No. Name Mixing Percentages of Components (558nm) SSA (558nm) 

5 sphericalmonabsorbing 1: 70%, 5: 30% 1 

51 spherical_med_dust 2: 72%, 6: 8%, 19: 20% 0.995 

70 spherical_med_coarse_dust 2: 20%, 19: 16%, 21: 64% 0.934 
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the mixing percentages of the M component aerosols involved. For each p, Xm=i 8pm = 1- 
MISR retrieval algorithm produces retrievals of AOD t p and compositions 6 P based on 
observed L p , independently for each pixel p on the spatial lattice, p = 1, . . . , P. 

Let T denote a pre-fixed grid of AOD values and s = 1, . . . , 21 denote the 21 MISR- 
designated single component aerosols. The Simulated MISR Ancillary Radiative Transfer 
(SMART) dataset contains the components {L RT ' s (r)}^l v r e T of radiation fields gener- 
ated by performing Radiative Transfer (RT) calculations on stratified atmospheric models 
containing each of the 21 component aerosols with AOD r [7]. Here we abuse the nota- 
tion by using M to denote the set of the M chosen component aerosols. Thus, for pixel p 
with a particular combination of its AOD value t p e T and compositions 6 P for M, define 
L rt (t p ,6 p ) as the model top-of-atmosphere sun reflectance: 

L RT (T p ,e p ) = Y J P sL RT ' s (T p ), 

seM 

a weighted linear combination of the radiation fields L RT ' S (-) obtained from SMART dataset. 
In order to retrieve AOD t p and compositions 6 P for pixel p, MISR retrieval algorithm 13 
first cross-generates a set of combinations of AOD values from a discrete grid and mixing 
percentages of component aerosols from the 74 MISR pre-fixed choices [6], namely 6 P e 
{0q, . . . , #q 4 }. Using as input each of the above cross-generated combinations of AOD and 
aerosol compositions, a Radiative-Transfer-Equation (RTE)-based simulation is conducted 
to calculate sun reflectance numerically. MISR's standard AOD retrieval algorithm is based 
on the comparison between such simulated radiances L RT and MISR observations L p = 
(Li p , ..., L Cp ) on C channels, with their goodness-of-fit measured by [6], 

X P = L 2^ ' (8) 

c=l L 



13 The model and algorithms in this study focus especially on aerosol retrievals over heterogenous land; 
retrieval over dark water is conceptually different and should be studied separately. 
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The combinations of AOD and compositions whose goodness-of-fit measures satisfy 
the pre-established threshold are considered successful retrievals and contribute to making 
inference of AOD values. 

7.3 Prior Distribution for AOD 

We model AOD's spatial dependence structure among adjacent pixels via a intrinsic GMRF 
prior. Define k as the homogenous scaler precision and Q > as the precision matrix. 
Assuming no bias in AOD's smoothness, r has its density in the form [23], 

MTk,e) = (2^)^ke|Jexp(-^T T er). (9) 

The precision matrix Q is symmetric and has the important property that, for pixels p' # p, 

Tp'-LTp | T- p ' p <==> Q p ' p = 0. 

In our model, for each pixel p, its AOD is directly related only to its few adjacent neighbor- 
hoods, suggesting a sparse precision matrix Q. Moreover, employing the intrinsic GMRF 
of first order [23], i.e. Q\ = 0, the density (9) is invariant to perturbation by the same 
constant to AOD of all pixels [3], 



P(t\k) oc k 2 exp 



Hs t '-4 



where ~ indicates adjacency on the two-dimensional lattice. In this case, the explicit form 
of precision matrix Q is, 

n p , if p'=p, 

Qp>p = '-\, if p'~ P , 

0, otherwise, 
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where n p is the number of adjacent pixels to pixel p;n p = 4 is all its neighboring pixels are 
cloudless. This allows us to model the AOD's smoothness regardless AOD's general level. 
The full conditional GMRF distribution of r is, 



_ ( n P K -AI2 



P(t p \t- p ,k) = (-^t) exp 



I Z n Pp':p'~P ) 



where r_ p are AOD values of all pixels other than the pth. 

7.4 Interpolation from SMART 

For each of the 21 component aerosols, denoted by s = 1, . . . , 21, MISR SMART dataset 
contains its RT-calculated radiances L RT ' S corresponding to AOD values on a fixed grid, 
i.e. (0, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1, 1.5, 2, 3, 4, 6). In order to obtain L RT ' S for an arbitrary 
AOD value between and 6, we consider L RT ' S as a function of AOD r and interpolate. 
Examples of L RT '\r) are shown in Figure 13, for Camera DF and blue band, Camera AN 
and blue band, Camera DF and red band respectively. An extensive exploratory data anal- 
ysis suggests a combination of linear and quadratic interpolations to balance between per- 
formance and computational efficiency. In particular, linear interpolations are used for 
r e (0, 0.4) U (3, 6), while quadratic interpolations are used for r 6 (0.4, 3). 

7.5 Sensitivity and Computational Expenses of Our MCMC Retrieval Algorithms 

We test the stability of outputs from the MCMC retrieval algorithm, using different initial 
values for the hyperparameters. As an example, results from two retrievals with initial val- 
ues for k = 100 and 200 are plotted in a scatterplot in the left panel of Figure 14. In the right 
panel plots retrievals with different initial or = (0.8,0.4,0.1,0.1) against (0.4,0.4,0.1,0.1). 
As is shown in Figure 14, the majority of the points scatter along the 45 degree line with 
few exceptions, indicating good stability behaviors of our algorithm outputs. 

The computation time of the three proposed algorithms are shown in Table 3, which 
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Figure 13: Examples of single- scattering radiances as a function of AOD, captured by 
different cameras and in different spectral bands. 
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Figure 14: Model sensitivity studies for k and a. 
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are all implemented on a server with 62.7 GB RAM using a typical MISR radiance im- 
agery. This case shows our MCMC algorithms in general improve the spatial resolution 
by a factor of 16. Though the parallel MCMC algorithm improves the MCMC algorithm's 
computational speed by a factor of four, it is still not efficient for online operation. There- 
fore, a Maximum a Posteriori (MAP) approach based on optimization is developed for 
operational purpose 14 . 



Algorithm 


Computation Time for Retrieval 


MISR (17.6 km) 


~ 16.7 minutes (1004.03 sec) 


Global MCMC (4.4 km) 


~ 6.2 hours (22452.26 sec) 


Parallel MCMC (4.4 km) 


- 1.6 hours (5863.77 sec) 


MAP (4.4 km) 


~ 6 minutes (378.91 sec) 



Table 3: Comparison of computation time. 



7.6 Specification of Posteriors and Proposal Distributions in the Metropolis-Hastings 



To sample t p (p = 1, . . . ,P) in the first P steps of the Gibbs sampler, we need to con- 
struct Markov chains with a stationary distribution as r^'s univariate marginal posterior 
distribution: 



| n„K 1 v-i ? v^ 
p(r p \r. p , 9, a, k, L) oc exp I — — (r p > iy ) - > 

I l "P p>:p>~p c=l 



(L cp - L c (t p , 6 P )) 
2^ 



We choose the following Gaussian proposal distribution, 



<5r(r*|r_ p ,/<-)ocexp 



1 p':p'~p J 



14 For the coherent structure of this study, the MAP retrieval algorithm is not included in this paper. 
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which is close to the marginal posterior distribution and easy to sample from [12]. After 
sampling r* from q(r*\T- p , k), the transition probability from iteration t to t + 1 is given by, 



Jf+i) _ 



t*, w.p. minll, 



Tp, otherwise. 



p(T*\T- p ,0,O;K,L)q(T { p'\T-p,K) 

p(jf \T- p ,0,<r,K,L)q(f\T- p ,k) 



When conditioning on a, r, and L, the compositions 6 P for each pixel p are independent 
and can be sampled separately. A Metropolis-Hastings algorithm is implemented with a 
Dirichlet proposal distribution: p(0* p \a) oc \\^ = i{6 pm ) am ~ l . 



1.1 Case Study for Improving Retrieval Resolution 

We take a closer look at a particular pixel and three of its adjacent retrievable pixels. Table 
4 shows the pixel for which our retrieval algorithm and MISR AOD retrievals disagree 
largely, with one of the largest errors. The center pixel is highlighted in gray, and is chosen 
due to the AERONET Beijing site within. Our Bayesian AOD retrieval is performed at 4.4 
km resolution; both retrievals are based on MISR observation over Beijing area (Location 

Table 4: Local comparison of MISR and our model 



Location 



1 

T 
T 



4 



AERONET 



0.2060 

N/A 
N/A 



N/A 



Our Retrievals MISR 



0.2168 
0.1870 
0.1389 



0.1540 
0.1540 
0.1540 



0.1710 



0.1540 



Path 122 Block 59) on January 26, 2009. While these four locations are represented by 
one pixel in MISR AOD algorithm at 17.6 km resolution and thus share the same retrieval 
output, our statistical retrieval algorithm at 4.4 km resolution produces four different AOD 
estimates respectively. 

Due to the existence of pixel 3, the special features of pixel 1 are averaged away with 
the relatively high AOD level, which could be one explanation as to why MISR's retrieval is 
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much lower than AERONET measurement in this particular example. This is one example 
of the benefits of including spatial smoothness into the model for a closer look at AOD 
spatial distribution. Certainly, such smoothness depends largely on the scale of observation 
and retrievals. Our algorithm will fail if certain events influence the local AOD level on 
a scale smaller than 4.4 km, which is not very uncommon. In that case, an algorithm that 
retrieves AOD at an even finer resolution is required to capture and model such impacts. 

7.8 More Details on Including Spatial Smoothness 



Bayesian AOD - MCMC 




Bayesian AOD - parallel without spatial smoothness 




I 



Figure 15: AOD retrievals with (top) and without (bottom) spatial smoothness. 

Introducing spatial smoothness into the retrieval model is seeking for a balance between 
particle movements in the open air and existence of local aerosol emissions. The exact 
calibration for such smoothness requires more studies on the speed and range of particulate 
matters' movements. However, Figure 15 comparing AOD retrievals with and without 
spatial smoothness, reveals the lack of an overall spatial structure leads to many cases 
of very sharp changes between adjacent pixels. Since such abrupt changes is not a very 
reasonable situation at this level of spatial resolution, the difference between the panels in 
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Figure 15 is possible results from lack of control for modeling and sampling errors, and 
confirms our argument for including spatial smoothness. In addition to the errors, retrieval 
results with a low level of smoothness are also confusing when trying to identify any major 
patterns of aerosol spatial distribution. 
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